df_apim <- df %>%# Reshaping to Actor-Partner format (4-field)reshape_dyadic_data(person_id ='userID',dyad_id ='coupleID',vars_to_reshape ='communication' ) %>%mutate(# Optional: grand-mean centeringcommunication_actor_gmc = communication_actor -mean(communication_actor),communication_partner_gmc = communication_partner -mean(communication_partner),# Create Dummy-Variables for male and femaleis_female =ifelse(gender ==1, 1, 0),is_male =ifelse(gender ==2, 1, 0), ) print_df(head(df_apim))
Distinguishable Dyads - Preparing Data
userID
coupleID
gender
is_male
is_female
communication
satisfaction
communication_actor
communication_partner
communication_actor_gmc
communication_partner_gmc
1_1
1
1
0
1
4.851107
4.841332
4.851107
6.004533
-0.1553604
0.9980657
1_2
1
2
1
0
6.004533
5.577165
6.004533
4.851107
0.9980657
-0.1553604
2_1
2
1
0
1
5.881321
7.248741
5.881321
6.433723
0.8748537
1.4272563
2_2
2
2
1
0
6.433723
6.960293
6.433723
5.881321
1.4272563
0.8748537
3_1
3
1
0
1
4.283971
6.928699
4.283971
5.516060
-0.7224960
0.5095933
3_2
3
2
1
0
5.516060
6.077688
5.516060
4.283971
0.5095933
-0.7224960
Distinguishable Dyads - Fitting the Model in BRMS
Model non-independence either via a shared random intercept or correlated residuals (not both). Posterior predictive checks, residual diagnostics and leave one out cross validation (e.g., via loo_compare()) may help inform which one to use.
For other families like ordinal regression with brms, one might work better than the other.
Distinguishable Dyads - Fitting the Model in BRMS
# Option A: Random-intercept model for non-independenceformula <-bf( satisfaction ~# Remove global intercept, introduce male and female intercepts.0+ is_male + is_female +# Actor effect for male and female communication_actor_gmc:is_male + communication_actor_gmc:is_female +# Partner effect for male and female communication_partner_gmc:is_male + communication_partner_gmc:is_female +# Option 0: Random intercept for male and female (correlated)# NOT identifiable in the cross-sectional case:# (0 + is_male + is_female | coupleID)# Option A: Shared dyad effect - interpretable dyad-level random intercept. (1| coupleID), # Two distinct residual variances for males vs females:sigma =~0+ is_male + is_female)priors <-c(# prior(normal(2, 3), class = "Intercept"),prior(normal(0, 5), class ="b"),prior(student_t(3, 0, 1.5), class ="b", dpar ="sigma"),prior(student_t(3, 0, 2.5), class ="sd"))model_dist_apim_a <-brm(formula = formula, data = df_apim,family =gaussian(link = identity),prior = priors,chains =4,cores =4,iter =2000,warmup =1000, seed =123,file =file.path('brms_cache', 'example1_dist_apim_a') # Cache the model)
Distinguishable Dyads - Fitting the Model in BRMS
# Option B: Residual CS for non-independence (no random intercept)formula <-bf( satisfaction ~# Remove global intercept, introduce male and female intercepts.0+ is_male + is_female +# Actor effect for male and female communication_actor_gmc:is_male + communication_actor_gmc:is_female +# Partner effect for male and female communication_partner_gmc:is_male + communication_partner_gmc:is_female +# Option B: Model non-independence by using compound symmetry of the residuals. # Aligns with the SEM conceptualization / literaturecosy(gr = coupleID),# Two distinct residual variances for males vs females:sigma =~0+ is_male + is_female)priors <-c(# prior(normal(2, 3), class = "Intercept"),prior(normal(0, 5), class ="b"),prior(student_t(3, 0, 1.5), class ="b", dpar ="sigma"),prior(beta(1, 3), class ="cosy"))model_dist_apim_b <-brm(formula = formula, data = df_apim,family =gaussian(link = identity),prior = priors,chains =4,cores =4,iter =2000,warmup =1000, seed =123,file =file.path('brms_cache', 'example1_dist_apim_b') # Cache the model)
Distinguishable Dyads - Check Model Convergence and Fit
Check Rhats and Effective Sample Sizes (ESS_tail and ESS_bulk) directly from the brms summary. Additionally you can (using Model B as an example):
0 of 4000 iterations saturated the maximum tree depth of 10.
Energy:
E-BFMI indicated no pathological behavior.
loo::pareto_k_table(loo(model_dist_apim_b))
All Pareto k estimates are good (k < 0.7).
Distinguishable Dyads - Check Model Convergence and Fit
plot(model_dist_apim_b, ask =FALSE)
Distinguishable Dyads - Check Model Convergence and Fit
pp_check(model_dist_apim_b, 'dens_overlay_grouped', group ='is_male')
Distinguishable Dyads - Check Model Convergence and Fit
pp_check(model_dist_apim_b, 'ecdf_overlay_grouped', group ='is_male')
Distinguishable Dyads - Check Model Convergence and Fit
pp_check(model_dist_apim_b, type ="loo_pit_overlay")
Distinguishable Dyads - Check Model Convergence and Fit
# Custom function to make DHARMa work with brms (see file 'Functions')DHARMa.check_brms(model_dist_apim_b)
Distinguishable Dyads - Check Model Convergence and Fit
In case of non-convergence, bad residuals or low predictive accuracy, the model is likely miss-specified.
You can try different families that match your data by using the family() argument. For instance, you can easily conduct an ordinal regression by simply setting family = cumulative() and adjusting your priors. For dichotomous outcomes you can use family = bernoulli(). Brms supports a variety of families and link functions for these families. In general, you should use the default link function (check ?brms::family.brmsfit()).
Note that for some of the later models in this tutorial, some ESS values are too low. In this case, usually we would increase iterations (and warmup), and/or adapt_delta (see brms documentation).
Distinguishable Dyads - Results
Est.
95% CI
Rhat
Bulk_ESS
Tail_ESS
Fixed Effects
is_male
4.60
[4.47, 4.74]
1.001
5163
3375
is_female
5.61
[5.50, 5.72]
1.000
4706
3152
is_male:communication_actor_gmc
1.84
[1.75, 1.94]
1.002
3380
3262
is_female:communication_actor_gmc
1.59
[1.51, 1.67]
1.002
3812
2935
is_male:communication_partner_gmc
0.17
[0.07, 0.27]
1.000
4128
3151
is_female:communication_partner_gmc
0.31
[0.23, 0.39]
1.000
3663
3475
Residual Structure
cosy
0.50
[0.44, 0.56]
1.001
4076
3169
sigma_is_male
0.47
[0.41, 0.53]
1.000
4362
2672
sigma_is_female
0.26
[0.20, 0.32]
1.001
4595
3241
Note that sigma is reported on the log scale (see ?brms::family.brmsfit()). Using exp(value) retrieves the simulated SD values well.
With repeated measures we need to center variables in order to not conflate levels of analysis:
Between-Person: The person-mean in relation to the grand mean
Within-Person: Daily fluctuations of an individual from their person-mean
While it may seem like we have 3-levels (couple/person/day), by including the means of both partners and their correlations in the model, all information about the couple-level (level 3) is included. This will become clear when comparing the exchangeable APIM to the DIM. We thus use a 2-level model with the dyad as the level of analysis.
Distinguishable Dyads - L-APIM - Centering
df_long_apim <- bmlm::isolate( df_long_apim, by ='userID', # NOT coupleIDvalue =c('provided_support_actor', 'provided_support_partner'),which ='both') %>%# renaming to avoid confusion with between- and within COUPLE centred DIM variablesrename(provided_support_actor_cwp = provided_support_actor_cw,provided_support_partner_cwp = provided_support_partner_cw,provided_support_actor_cbp = provided_support_actor_cb,provided_support_partner_cbp = provided_support_partner_cb, ) %>%# Centering timemutate(diaryday_c =scale(diaryday, center =TRUE, scale =FALSE) )
Distinguishable Dyads - L-APIM - Model nlme
nlme model of this data for reference (simplified random effects structure to achieve convergence).
library(nlme)model1 <-lme(fixed =# Intercepts closeness ~0+ is_male + is_female +# Time Slopes is_male:diaryday_c + is_female:diaryday_c +# Between-Person APIM is_male:provided_support_actor_cbp + is_male:provided_support_partner_cbp + is_female:provided_support_actor_cbp + is_female:provided_support_partner_cbp +# Within-Person APIM is_male:provided_support_actor_cwp + is_male:provided_support_partner_cwp + is_female:provided_support_actor_cwp + is_female:provided_support_partner_cwp,random =~0+ is_male + is_female + is_male:diaryday_c + is_female:diaryday_c# in nlme we need to remove many random slopes to achieve convergence# alternatively, we could try to remove correlations between the slopes. #+ is_male:provided_support_actor_cwp #+ is_male:provided_support_partner_cwp#+ is_female:provided_support_actor_cwp #+ is_female:provided_support_partner_cwp | coupleID, weights =varIdent(form =~1| gender), # heterogeneous residual variancescorr =corCompSymm(form =~1| coupleID/diaryday), # compound symmetry (partner's values on each day)data = df_long_apim, control =list(maxIter=1000))summary(model1)
Distinguishable Dyads - L-APIM - Model brms
Equivalent brms model with simplified random effects structure.
formula <-bf( closeness ~0+ is_male + is_female +# Time Slopes is_male:diaryday_c + is_female:diaryday_c +# Between-Person APIM is_male:provided_support_actor_cbp + is_male:provided_support_partner_cbp + is_female:provided_support_actor_cbp + is_female:provided_support_partner_cbp +# Within-Person APIM is_male:provided_support_actor_cwp + is_male:provided_support_partner_cwp + is_female:provided_support_actor_cwp + is_female:provided_support_partner_cwp +# Accounting for non-independence between partners' means and trajectories # and effect sensitivities via random effects: (0+ is_male + is_female+ is_male:diaryday_c + is_female:diaryday_c| coupleID) +# Accounting for daily non-independence (mimicking nlme's corCompSym)# modelled via common 'shocks' with a random intercept: (1| coupleID:diaryday)# heteroscedastic residuals , sigma ~0+ is_male + is_female )priors <-c(prior(normal(4, 2), class ="b", coef ="is_male"), # male interceptprior(normal(4, 2), class ="b", coef ="is_female"), # female interceptprior(normal(0, 2), class ="b") # other beta coefficients# Other priors can be specified, but brms choosed defaults well. # Defaults can be seen via default_prior(formula, data))model_dist_apim_long_simple <-brm(formula = formula, data = df_long_apim,family =gaussian(link = identity), # student() is also often a nice robust option. prior = priors,chains =4,cores =4,iter =2000,warmup =1000, seed =123,file =file.path('brms_cache', 'example1_dist_apim_long_simpel') # Cache the model)summary(model_dist_apim_long_simple)
... +ar(time = diaryday, gr = coupleID:userID, p =1), sigma ~0+ is_male + is_female)
Distinguishable Dyads - L-APIM - AR1 Model
Warning: Often, the daily shocks and the time slopes are enough and adding AR1 additionally may “soak up” too much variance. this can bias fixed effects estimates towards 0.
You may want to check which model performs better!
The difference in predictive accuracy, as indexed by Expected Log Predictive
Density (ELPD-LOO), suggests that 'model_dist_apim_long_simple' is the best
model (ELPD = -18689.18), followed by 'model_dist_apim_long_ar' (diff-ELPD =
-107.87 +- 25.27, p < .001)
In this instance, the AR1 model performs much worse, indicating potential overfitting or other issues. Even in the case of no clear difference, it may be advisable to choose the more parsimonious model, which is the one without AR1.
Maximal Model Specification (no AR(1))
We add all random slopes and their correlations. This is not converging in nlme.
formula <-bf( closeness ~# Intercepts0+ is_male + is_female +# Time Slopes is_male:diaryday_c + is_female:diaryday_c +# Between-Person APIM is_male:provided_support_actor_cbp + is_male:provided_support_partner_cbp + is_female:provided_support_actor_cbp + is_female:provided_support_partner_cbp +# Within-Person APIM is_male:provided_support_actor_cwp + is_male:provided_support_partner_cwp + is_female:provided_support_actor_cwp + is_female:provided_support_partner_cwp +# Accounting for non-independence between partners' means and trajectories # and effect sensitivities via random effects: (0+ is_male + is_female + is_male:diaryday_c + is_female:diaryday_c + is_male:provided_support_actor_cwp + is_male:provided_support_partner_cwp + is_female:provided_support_actor_cwp + is_female:provided_support_partner_cwp | coupleID ) +# Accounting for daily non-independence (1| coupleID:diaryday), # No more AR1. You could also test MA or ARMA. sigma ~0+ is_male + is_female # heteroscedastic residuals)priors <-c(prior(normal(4, 2), class ="b", coef ="is_male"),prior(normal(4, 2), class ="b", coef ="is_female"),prior(normal(0, 2), class ="b"))model_dist_apim_long_complex <-brm(formula = formula, data = df_long_apim,family =gaussian(link = identity),prior = priors,chains =4,cores =4,iter =2000,warmup =1000, seed =123,file =file.path('brms_cache', 'example1_dist_apim_long_complex') # Cache the model)
In complex models like this, we want to make sure we are not overfitting.
Conduct regular model convergence checks (inspecting chains, Rhats, ESS, Multicollinearity, Residual diagnostics etc.). If these are bad (but good for a simpler model), it could be a sign of overfitting.
Leave-one-out cross-validation is a powerful tool for investigating overfitting. It tests out of sample prediction and thus punishes overfitting.
In case of High Pareto-Ks there are influential datapoints present. This can be a sign of overfitting but it is not necessarily a problem. But it is recommended to use loo(model, moment_match = TRUE) in these cases (the function will tell you).
The loo values (elpd) themselves cannot be interpreted individually, but we can compare models again (next page)
No obvious overfitting detected. In fact, the more complex model seems to be more predictive. Some rule of thumbs say that if the difference is higher than 2x the SE or 4x the SE we can say there is a difference.
If there was no obvious difference, we could follow our theory and the philosophy of the maximal random effects structure (Barr et al., 2013). Alternatively, we could follow the philosophy of parsimony, especially in the territory of such complex models.
(continued on next page)
Maximal Model Specification: Diagnostics
One approach is to use the report package’s approach to obtain a p-value for the comparison (see their documentation for which assumptions they make):
report::report(a)
The difference in predictive accuracy, as indexed by Expected Log Predictive
Density (ELPD-LOO), suggests that 'model_dist_apim_long_simple' is the best
model (ELPD = -18689.18), followed by 'model_dist_apim_long_complex' (diff-ELPD
= -3.87 +- 1.59, p = 0.015)
In this instance, the maximal model (no AR) seems to be the preferred one.
In the case of overfitting, we could try to improve the model by subsequently excluding the smallest random effects (or correlations between them) and compare various models. For more guidance on random effects check out del Rosario & West (2025).
Maximal Model Specification: Comparing Estimates to Simulated Values (Ground Truth)
Same cross-sectional data in the same 4-field actor-partner format as before.
print_df(head(df_apim))
userID
coupleID
gender
is_male
is_female
communication
satisfaction
communication_actor
communication_partner
communication_actor_gmc
communication_partner_gmc
1_1
1
1
0
1
4.851107
4.841332
4.851107
6.004533
-0.1553604
0.9980657
1_2
1
2
1
0
6.004533
5.577165
6.004533
4.851107
0.9980657
-0.1553604
2_1
2
1
0
1
5.881321
7.248741
5.881321
6.433723
0.8748537
1.4272563
2_2
2
2
1
0
6.433723
6.960293
6.433723
5.881321
1.4272563
0.8748537
3_1
3
1
0
1
4.283971
6.928699
4.283971
5.516060
-0.7224960
0.5095933
3_2
3
2
1
0
5.516060
6.077688
5.516060
4.283971
0.5095933
-0.7224960
Exchangeable Dyads - Cross-Sectional APIM: Model
formula <-bf( satisfaction ~1+ communication_actor_gmc + communication_partner_gmc +# Option 1: A single couple level random intercept.# In the cross sectional case this is the maximal random effects structure and# (1 | coupleID)# Option 2: using a compound symmetry residual structure. cosy(gr = coupleID)# Note: no need to model separate sigmas for each partner.# Homogeneous residual variance is estimated. # Implied: sigma = ~ 1)priors <-c(prior(normal(2, 10), class ="Intercept"),prior(normal(0, 5), class ="b"),prior(student_t(3, 0, 1.5), class ="sigma"),prior(beta(1, 3), class ="cosy"))model_ind_apim <-brm(formula = formula, data = df_apim,family =gaussian(link = identity),prior = priors,chains =4,cores =4,iter =2000,warmup =1000, seed =123,file =file.path('brms_cache', 'example1_ind_apim') # Cache the model)
Exchangeable Dyads - Cross-Sectional APIM: Test for Distinguishability
Leave-one-out (loo) cross-validation for model comparison (even if not nested).
a <-loo_compare(loo(model_ind_apim), loo(model_dist_apim_b))print(a)
Model
elpd_diff
se_diff
model_dist_apim
0.0
0.0
model_ind_apim
-172.9
16.5
report::report(a)
The difference in predictive accuracy, as indexed by Expected Log Predictive Density (ELPD-LOO), suggests that ‘model_dist_apim_b’ is the best model (ELPD = -1806.42), followed by ‘model_ind_apim’ (diff-ELPD = -172.85 +- 16.53, p < .001). See: documentation of “report”
Top sliders: grand-mean-centered Communication for Actor and Partner (APIM).
Bottom sliders: reparameterized communication — Centered Between-Couple (xcbc) and Centered Within-Couple (xcwc) (DIM).
Equivalence APIM and DIM
If the couple mean goes up by 1 and the within-couple deviation from the couple mean stays fixed, this means that both partners’ predictors must go up by 1 unit. Thus, the effects are linear combinations:
\[b_{cbc} = b_{actor\_gmc} + b_{partner\_gmc}\]
Similarly, if a person’s deviation from their couple mean is one unit higher and the couple mean stays fixed, this means their partners’ deviation must be one unit lower. Thus, the effects are a linear combination:
Exchangeable Dyads - Longitudinal APIM/DIM: Simulated Data
We use the same data as in the exchangeable case. But we don’t need the gender-related variables anymore, as we assume that partners are indistinguishable.
userID
coupleID
diaryday
diaryday_c
closeness
provided_support_actor
provided_support_partner
provided_support_actor_cbp
provided_support_actor_cwp
provided_support_partner_cbp
provided_support_partner_cwp
31_1
31
0
-27
5.47
6.02
7.09
0.65
0.41
0.58
1.54
31_1
31
1
-26
7.69
6.13
6.95
0.65
0.51
0.58
1.4
31_1
31
2
-25
5.83
6.1
6.34
0.65
0.49
0.58
0.79
31_1
31
3
-24
6.55
6.99
5.39
0.65
1.38
0.58
-0.15
...
...
...
...
...
...
...
...
...
...
...
31_2
31
0
-27
5.19
7.09
6.02
0.58
1.54
0.65
0.41
31_2
31
1
-26
6.96
6.95
6.13
0.58
1.4
0.65
0.51
31_2
31
2
-25
6.9
6.34
6.1
0.58
0.79
0.65
0.49
31_2
31
3
-24
7.59
5.39
6.99
0.58
-0.15
0.65
1.38
Exchangeable Dyads - Longitudinal DIM/APIM
We need an APIM or DIM on both the between person level and the within-person level.
Due to equivalence, we could have a DIM on one level and an APIM on the other.
If we want a within-person DIM:
df_long_dim <- df_long_apim2 %>%mutate(# Between Person Level DIM (centering between couples and within couple between person)provided_support_cbc = (provided_support_actor_cbp + provided_support_partner_cbp) /2,provided_support_cwcbp = (provided_support_actor_cbp - provided_support_partner_cbp) /2,# Within Person Level DIM (couple mean and deviation on each day)provided_support_cwp_mean = (provided_support_actor_cwp + provided_support_partner_cwp) /2,provided_support_cwp_halfdiff = (provided_support_actor_cwp - provided_support_partner_cwp) /2 )
Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the sum and difference approach
We need to randomly assign -1 and +1 for each Partner within each couple. This will be needed for contrast coding the intercept.
For example:
df_long_dim <- df_long_dim %>%group_by(coupleID) %>%mutate(base =ifelse(userID ==min(userID), 1, -1),flip =1-2*rbinom(1, 1, 0.5), # yields +1 or -1 once per coupleIdiff = base * flip ) %>%ungroup() %>%relocate(Idiff, .after = coupleID) %>% dplyr::select(-base, -flip)print_df(print_couple_preview(df_long_dim, couple_id ="31"))
Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the Sum and Difference Approach
userID
coupleID
Idiff
diaryday
diaryday_c
closeness
provided_support_actor
provided_support_partner
provided_support_actor_cbp
provided_support_actor_cwp
provided_support_partner_cbp
provided_support_partner_cwp
provided_support_cbc
provided_support_cwcbp
provided_support_cwp_mean
provided_support_cwp_halfdiff
31_1
31
-1
0
-27
5.47
6.02
7.09
0.65
0.41
0.58
1.54
0.61
0.04
0.97
-0.57
31_1
31
-1
1
-26
7.69
6.13
6.95
0.65
0.51
0.58
1.4
0.61
0.04
0.96
-0.44
31_1
31
-1
2
-25
5.83
6.1
6.34
0.65
0.49
0.58
0.79
0.61
0.04
0.64
-0.15
31_1
31
-1
3
-24
6.55
6.99
5.39
0.65
1.38
0.58
-0.15
0.61
0.04
0.61
0.77
31_1
31
-1
4
-23
7.03
6.7
5.19
0.65
1.08
0.58
-0.36
0.61
0.04
0.36
0.72
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
31_2
31
1
0
-27
5.19
7.09
6.02
0.58
1.54
0.65
0.41
0.61
-0.04
0.97
0.57
31_2
31
1
1
-26
6.96
6.95
6.13
0.58
1.4
0.65
0.51
0.61
-0.04
0.96
0.44
31_2
31
1
2
-25
6.9
6.34
6.1
0.58
0.79
0.65
0.49
0.61
-0.04
0.64
0.15
31_2
31
1
3
-24
7.59
5.39
6.99
0.58
-0.15
0.65
1.38
0.61
-0.04
0.61
-0.77
31_2
31
1
4
-23
7.55
5.19
6.7
0.58
-0.36
0.65
1.08
0.61
-0.04
0.36
-0.72
Exchangeable Dyads - Longitudinal APIM/DIM: Model
formula <-bf( closeness ~1+ diaryday_c +# Within-person APIM provided_support_actor_cwp + provided_support_partner_cwp +# Equivalent to within-person DIM:# provided_support_cwp_mean + # provided_support_cwp_halfdiff +# Between-person APIM provided_support_actor_cbp + provided_support_partner_cbp +# Equivalent to between-person DIM# provided_support_cwcbp+# provided_support_cbc +# Dyad-Level intercept and slopes for time-varying predictors (1+ diaryday_c + provided_support_actor_cwp + provided_support_partner_cwp | coupleID) +# Both partners' deviations from these dyad-level means and slopes.# Note: separate random effects block to make them uncorrelated form dyad-level REs. (0+ Idiff +I(Idiff * diaryday_c) +I(Idiff * provided_support_actor_cwp) +I(Idiff * provided_support_partner_cwp) | coupleID) +# Accounting for same-day shocks/coupling (1| coupleID:diaryday)# Autocorrelated residuals# Due to bad fit in the distinguishable case, we omit AR1 from the start.# ar(time = diaryday, gr = coupleID:userID, p = 1)# Again, no need to model heterogeneous residual variances (sigma)# Implied: sigma = ~ 1)priors <-c(prior(normal(4, 2), class ="Intercept"),prior(normal(0, 2), class ="b")# We could set priors for other things, but brms sets good default priors#prior(student_t(3, 0, 1.5), class = "sd"),#prior(lkj(2), class = "cor"), # correlation prior for random effect matrix#prior(student_t(3, 0, 1.5), class = "sigma"))model_apim_ind_long <-brm(formula = formula, data = df_long_dim,family =gaussian(link = identity),prior = priors,chains =4,cores =4,iter =2000,warmup =1000, seed =123,file =file.path('brms_cache', 'model_apim_ind_long_apim'))
Randomly give one partner a constant of -1 and the other a constant of 1
The couple level intercept represents the mean (or sum) of both partners’ intercepts. (1 | coupleID)
A column of 1s and -1s represent deviations (difference) of each partner from the couple intercept with each partner contributing equally in one direction (+1 and -1)
Exchangeability is retained: if partners are flipped, the result is the same.
(del Rosario & West (2025) provide practical guidance on how to reduce the random effects structure in case of non-convergence.)
Extract Full APIM Random Effects Variance-Covariance Matrix
We can rotate the random effects structure back to obtain the full APIM Random effects variance-covariance matrix (Kenny & Ackerman, 2023). Just as we would in a SEM model with equality constraints!
For example, the within-person variance of the intercept is:
\[var(Intercept) + var(Idiff)\]
and the cross-person covariance of both partners’ intercepts is:
\[var(Intercept) - var(Idiff)\]
Note: This only holds true when the means and Idiffs are uncorrelated. This is why we put them in separate random effects blocks.
Extract Full APIM Random Effects Variance-Covariance Matrix
Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278. https://doi.org/10.1016/j.jml.2012.11.001
Bolger, N., Laurenceau, J.-P., & DiGiovanni, A. (2025). Unified analysis model for indistinguishable and distinguishable dyads. Innovations in Interpersonal Relationships and Health Research: Advancing the Integration of Interdisciplinary Approaches to Dyadic Behavior Change. https://doi.org/10.17605/OSF.IO/WYDCJ
del Rosario, K. S., & West, T. V. (2025). A Practical Guide to Specifying Random Effects in Longitudinal Dyadic Multilevel Modeling. Advances in Methods and Practices in Psychological Science, 8(3), 25152459251351286. https://doi.org/10.1177/25152459251351286
Gabry, J., Češnovar, R., Johnson, A., & Bronder, S. (2024). Cmdstanr: R interface to CmdStan. https://mc-stan.org/cmdstanr/
Guo, J., Gabry, J., Goodrich, B., Johnson, A., Weber, S., & Badr, H. S. (2025). Rstan: R interface to stan. https://mc-stan.org/rstan/
Iida, M., Seidman, G., & Shrout, P. E. (2018). Models of interdependent individuals versus dyadic processes in relationship research. Journal of Social and Personal Relationships, 35(1), 59–88. https://doi.org/10.1177/0265407517725407
Kenny, D. A., & Ackerman, R. A. (2023). Estimation of random effects with over-time dyadic data using multilevel modeling: The sum and difference method. OSF Preprints. https://doi.org/10.31219/osf.io/fju72
Kenny, D. A., & Cook, W. (1999). Partner effects in relationship research: Conceptual issues, analytic difficulties, and illustrations. Personal Relationships, 6(4), 433–448. https://doi.org/10.1111/j.1475-6811.1999.tb00202.x
Kenny, D. A., & Kashy, D. A. (2011). Dyadic Data Analysis Using Multilevel Modeling. In Handbook of Advanced Multilevel Analysis. Routledge.
Kenny, D. A., Kashy, D. A., & Cook, W. L. (2006). Dyadic data analysis. Guilford Press.
Lüdecke, D., Makowski, D., Ben-Shachar, M. S., Patil, I., Wiernik, B. M., Bacher, E., & Thériault, R. (2025). Easystats: Framework for easy statistical modeling, visualization, and reporting. https://easystats.github.io/easystats/
Stadler, G., Scholz, U., Bolger, N., Shrout, P. E., Knoll, N., & Lüscher, J. (2023). How is companionship related to romantic partners’ affect, relationship satisfaction, and health behavior? Using a longitudinal dyadic score model to understand daily and couple-level effects of a dyadic predictor. Applied Psychology: Health and Well-Being, 15(4), 1530–1554. https://doi.org/10.1111/aphw.12450
Vehtari, A., Gabry, J., Magnusson, M., Yao, Y., Bürkner, P.-C., Paananen, T., & Gelman, A. (2024). Loo: Efficient leave-one-out cross-validation and WAIC for bayesian models. https://mc-stan.org/loo/